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The  models  proposed  should  prove  to  be  useful  for 
planning  and  evaluating  both  new  and  existing  time-sharing 


computer  systems. 
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2.  Model  1:  Diffusion  Approximations  for  One  Job  Type 


FIGURE  1 


In  this  section  we  assume  that  each  of  N terminals  is 
active  auid  requests  service  according  to  an  exponential  distribution 
with  parameter  X . Each  job  requires  service  from  the  computer 
and  the  service  times  also  have  an  exponential  distribution  with 
parameter  Nu#  independently  of  the  arrivals.  All  jobs  are 
identical,  and  the  computer  is  assumed  to  operate  according  to 
a first-coroe,  first-served  (FCFS)  queue  discipline. 

We  define  Q{t)  to  be  the  number  of  jobs  at  the  computer 
at  time  t.  Clearly  {Q(t) , t ^ 0}  is  a continuous-time  Markov 
chain.  One  can  explicitly  solve  the  steady  state  birth-death 
equations  (see  Karlin  (1967,  p.208)  for  example).  This  model 
has  been  extensively  studied  as  the  "repairman"  model;  Feller 
(1957)  . 


It  is  somewhat  difficult  to  understand  the  dynamic  or 
transient  behavior  of  such  a system.  One  could  solve  the 
Kolmogorov  forward  equations,  perhaps  in  closed  form  by  means  of 
Laplace  transforms,  or  perhaps  numerically;  however,  an  approach 
which  gives  a great  deal  of  insight  into  the  dynaunics  of  the 
system  is  to  introduce  a diffusion  approximation  for  {Q(t) , t > 0} 
We  note  that  if  the  number  of  terminals  N is  large,  changes 
in  Q(t)  will  occur  frequently.  Over  any  short  period  of  time, 
the  system  size  will  change  by  the  sum  of  many  independent 
Bernoulli-like  random  variables.  In  the  spirit  of  the  central 
limit  theorem,  then,  we  model  these  changes  by  normal  random 
varied>les.  The  resulting  state  space  will  be  ]R  ^ = [0,“) 
rather  than  {0,  1,  ...  , N).  A rigorous  theory  of  this  approxi- 
mation has  been  supplied  by  Barbour  (1974) • 

The  process  {Q(t) , t ^ 0}  is  thus  approximated  by  a 
continuous-time  Markov  process  satisfying  the  Ito  type  stochastic 
differential  equation 

dQ(t)  - X(N  - Q(t)]dt  - pN  dt 

+ /X(N  - Q(t))  dWj^(t)  - /wN  dW^(t)  (2.1) 

where  W^(t)  and  W^(t)  are  independent  standard  Wiener 
processes.  Actually,  the  parameters  X and  u may  be  time 
dependent.  For  what  follows  we  assume  X and  p constant, 
and  y/X  < 1.  The  argument  for  the  model  (2.1)  is  simply  as 
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follows.  The  term  X [N  - Q(t)]  represents  the  expected  rate  of 
arrivals  when  the  number  of  users  at  the  terminal  is  Q(t) , and 
uN  represents  the  corresponding  expected  rate  of  departures; 
together  they  are  called  the  drift.  The  term  /X  [N-Q (t) ) dWj^(t) 
represents  the  random  flucutation  of  the  arrival  component  to 
Q*s  increase,  and  /yN  dW2(t)  represents  the  corresponding 
independent  fluctuation  in  departures.  Our  approximation  regards 
both  of  these  fluctuations  as  independently  Gaussian. 

Next  we  define  a new  process  {Xjj(t),  t ^ 0}: 

X^(t)  = (Q(t)  -Nq(t))//N  or  Q(t)=Nq(t)  +/NXjj(t)  . (2.2) 

Using  Ito's  Lemma  (Arnold,  p.  90),  one  may  derive  the 
stochastic  differential  equation  governing  (X^^Ct),  t ^ 0}. 

We  find 

dXj^{t)  - /N  ^ ^1^-'  - X(l-q{t))jdt  + y dt  - XXj^{t)dt 

+ ✓X(l-q(t))  dWj^(t)  + dW2(t)  (2.3) 

Now  let  N -►  + ».  For  {Xjj(t),  t ^ 0}  to  converge  weakly 
to  a finite  limit  {X(t),  t > 0} , the  coefficient  of  the  /n  dt 
term  must  be  identically  zero  for  all  t.  This  compels  us  to 
state  that 
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q'(t)  = A(l-q(t)) 


9 


(2.4) 


- U 


q(0) 


= Q(0) 

N 


Equation  (2.4)  is  easily  solved  to  give 

q(t)  = (q(0)  - (l-w/X))e’^^  + (1-u/A),  M < X (2.5) 

When  li  ^ X this  particular  approximation  has  little  value,  for 
it  predicts  an  eventual  queue  of  zero  length;  this  occurs  because 
we  have  not  accounted  for  the  presence  of  a boundary  at  Q = 0, 
encountered  with  noticeable  frequency  only  when  the  system  is  very 
lightly  loaded.  Note  that  in  such  a situation  the  appropriate 
diffusion  approximation  is  that  for  a single-server  queue  with 
Poisson  arrivals;  see  Gaver  (1968) . Here  we  model  systems  that 
are  rather  heavily  loaded,  and  ignore  the  boundary. 

It  q(t)  is  chosen  as  in  (2.5),  then  the  results  of 
Kurtz  (1971)  and  Barbour  (1974)  insure  that  {Xjj(t),  t ^ 0} 
will  converge  wea)cly  to  {X(t),  t ^ 0},  a diffusion  process 
governed  by  the  stochastic  differential  equation 

dX(t)«-XX(t)  dt  + /X(l-q(t)  ) + v dW(t)  , X(0)  = 0.  (2.6) 

Note  that  we  have  consolidated  the  two  random  fluctuation  terms 
into  one. 

Equation  (2.6)  characterizes  a non-stationary  Ornstein- 
Uhlenbec)c  process.  Many  results  are  )cnown  for  these  processes 
(see  Arnold,  Chapter  8).  The  diffusion  approximation  of 
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{Q(t) , t ^ 0}  amounts  to  approximating  Q(t)  by 
Nq(t)  + X(t)  rather  than  Nq(t)  + Xjj(t).  The  quantity 
Nq(t)  is  referred  to  as  the  deterministic  approximation, 
while  i/Tj  X(t)  is  a stochastic  noise  process  superimposed  on 
the  deterministic  term. 

One  can  use  (2.5)  and  (2.6)  to  study  the  dynamic  or 
transient  behavior  of  the  system.  If  we  let  t an  account 

of  the  steady-state  behavior  is  obtained.  In  all  cases  Q(t) 
will  be  approximately  normally  distributed.  The  parameters  of 
this  distribution  are  easily  characterized.  If  v(t)  = Var(X(t)), 
then  (Arnold,  Chapter  8) 


E(X(t))  = 0 

and 

V' (t)  = -2X(t)  + A(l-q(t))  + u,  v(0)  = 0 
Equation  (2.7)  can  be  easily  solved  using  (2.5)  to  give 

v(t)  = ^ (1  - e"^^)  + (1  - q(0))(e’^^  - e'^^^)  . (2.8) 

A 

The  marginal  distribution  of  Q(t)  is  approximately 
normal  with  mean  Nq(t)  and  variance  Nv(t) . The  transient 
behavior  can  be  determined  using  the  above  analysis.  If  at  some 
time  tg , Q(tQ)  * Nq^,  then  Q(t)  will  be  approximately  normally 
distributed  with  mean  and  variance  given  by  (2.5)  and  (2.8)  with 
t replaced  by  ^"^0  q(0)  replaced  by  q^ . 
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As  t , steady  state  prevails.  The  distribution  of 

Q(t)  converges  weakly  to  a normal  distribution  with  mean 
N(1  - \i/\)  and  variance  Ny/A. 

Equivalent  diffusion  approximation  ideas,  but  given  in 
a much  different  form,  were  originally  presented  by  Iglehart 
(1965)  . 
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3 . Model  2;  Two  Types  of  Customers 


Model  I is  clearly  an  oversimplified  representation  of 
a time-sharing  computer  system.  The  computer,  considered  as  a 
server,  actually  performs  complex  functions,  and  appears  to  be 
a network  of  queues  and  servers.  We  will,  throughout  this 
paper,  continue  to  model  the  computer  as  a single  server  with 
a f irst-come-f irst-served  queue  discipline.  However,  we  will 
generalize  the  exponential  service  distribution  to  a wider 
class,  say  the  class  of  phase  distributions.  A second  problem 
arises  with  the  assumption  of  exponential  user  "think  times," 
which  again  may  appear  unrealistic.  Finally,  it  is  unreason- 
able to  assume  that  all  of  the  terminal  users  request  the 
same  type  of  service  from  the  computer.  In  fact,  some  users 
may  be  performing  simple  editing  tasks,  others  may  be  making 
heavy  use  of  the  computer  for  scientific  computations,  while 
others  may  be  using  it  in  a mixture  of  these  ways.  This  type 
of  situation  with  heterogeneous  job  types  is  often  modelled 
by  assuming  that  the  effective  service  time  is  a fixed  mixture 
of  individual  service  distributions.  Unfortunately,  in  this 
context  such  an  assumption  seems  unsatisfactory,  since  the 
appropriate  mixture  must  ultimately  depend  on  the  types  of 
jobs  in  service  and  must  therefore  change  dynamically. 

In  modelling  computer  systems  with  a mixture  of  job 
types,  a processor  sharing  queue  discipline  is  typically 
assumed.  Processor  sharing  is  the  limiting  version  of  a 
"round  robin"  queue  discipline  as  the  quantum  of  service  given 


10 


to  each  job  approaches  zero.  Processor  sharing  provides 
great  mathematical  cc  venience  in  that  resulting  queueing 
networks  will,  unde  veak  assumptions,  exhibit  a product-form 
steady  state  distribution  (see  Baskett  et  al.  (1975)). 
Unfortunately,  processor  sharing  may  not  faithfully  model  the 
behavior  of  the  computer.  The  computer  scheduling  policy 
is  closer  to  a FCFS  queueing  system  than  to  a processor 
sharing  one,  assuming  that  actual  quanta  are  relatively 
large,  and  most  jobs  finish  service  before  the  end  of  a 
quantum.  The  models  that  we  consider  better  portray  a batch 
system  with  terminal  users  than  a true  time-sharing  system. 
The  central-server  type  of  multiprogramming  model  will  be 
considered  in  another  paper. 

In  this  section,  we  assume  there  are  two  types  of 
jobs,  each  with  an  exponential  service  time,  with  parameters 
and  ^2*  Model  2 we  assume  that  each  time  a job  is 

submitted  from  a terminal,  there  is  a probability  p of  its 
being  type  l,and  1-p  of  its  being  Type  2.  We  assume  jobs 
arrive  from  any  terminal  independently  according  to  an 
exponential  (A)  distribution  and  that  the  computer  handles 
jobs  in  a FCFS  fashion. 

This  system  can  be  simply  modelled  as  a Markov  chain. 
Ordinarily  with  two  types  of  customers  in  a single  server 
queue,  one  must  keep  track  of  the  order  of  the  jobs  in  the 
queue.  In  this  case  , however,  customers  need  not  have  their 
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identities  established  until  they  enter  service.  Consequently, 
if  Q(t)  represents  the  number  of  jobs  in  service  at  time  t 
and  I(t)  = i if  a type  i job  is  in  service,  then 
{ (N (t) , I (t) ) , t ^ 0}  is  a continuous  time  Markov  chain  with 
transitions  given  by 


Transition 

(M,l ) — > (M+1,1)  M > 0 

(M-1,1)  M > 1 

(M-1,2)  M > 1 

(M,2)  — > (M+1,2)  M > 0 

(M-1,1)  M ^ 1 

(M-1,2)  M > 1 

(0,1)  — > (1,1) 

(0,2)  — > (1,2) 


Rate 

X(N-M)  dt  + o(dt) 

Pj^p  dt  + o(dt) 

(1-p)  dt  + o (dt) 

A(N-M)  dt  + o(dt) 
y^p  dt  + o(dt) 

U2 (1-p)  dt  + o (dt) 

(3.1) 

XNp  dt  + o(dt) 

XN(l-p)  dt  + o(dt) 


The  states  (0,1)  and  (0,2)  can  be  aggregated  into  a 
single  state  representing  system  emptiness. 

The  steady  state  equations  can  be  solved,  although  not 
in  a convenient  closed  form.  A Gauss-Seidel  iteration  procedure 
will  provide  an  efficient  numerical  approach  even  for  very 
large  N in  view  of  the  sparseness  of  the  transition  matrix. 
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We  wish  to  study  three  approximation  methods  for  this 
simple  model.  Based  on  their  performance  we  may  be  able  to 
use  them  in  more  complicated  models.  The  first  method  is  a 
straightforward  generalization  of  Model  1 and  is  based  heavily 
on  the  fact  that  any  terminal  can  submit  either  type  of  job; 
an  assumption  requiring  jobs  to  be  typed  only  upon  reaching 
service.  To  approximate  the  behavior  of  this  system  in  the 
multitype  setting,  we  can  replace  the  two  service  rates  p 
and  ^2  ® weighted  rate  p*  = (p/y^^  + (l-p)/P2)  This 

approximation  is  intuitive  in  that  it  gives  the  correct  mean 
service  rate  and  is  easy  to  handle.  One  can  apply  the  analysis 
given  in  Section  2 by  merely  replacing  p by  p*.  Numerical 
examples  of  the  accuracy  of  this  method  will  be  .given  later. 

In  general,  it  provides  a high  level  of  accuracy  when  compared 
with  an  exact  solution  using  (3.1),  but  the  reader  should 
note  that  it  will  not  be  easily  applied  to  situations  in 

which  the  number  of  terminals  of  each  type  is  fixed  and  each 
terminal  submits  jobs  of  a particular  type  only. 

The  second  approximation  is  based  upon  a processor 
sharing  queue  discipline  concept  modified  to  represent  FCFS 
system  behavior.  If  the  queue  were  to  contain  N^(t)  jobs  of 
type  i at  time  i,  i = 1,  2,  then  the  traditional  processor 
sharing  model  allocates  N^(t)/{N^(t)  + N2(t))  time  units  to 
jobs  of  type  i,  i = 1,  2 during  each  (t,  t + dt)  time 
interval.  Each  job  waiting  for  service  is  given  equal  weight; 
consequently,  short  jobs  are  given  preferential  treatment. 
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for  they  are  not  blocked  behind  long  jobs.  For  this  queue 
discipline  {(Nj^(t),  N2(t)),  t > 0}  is  a continuous  time 
Markov  chain  and,  using  the  results  of  Baskett  et  al.  (1975), 
the  stationary  distribution  will  be  of  a product  form.  A 
simple  modification  will  allow  us  to  model  the  behavior  of 
a FCFS  server.  We  weight  each  job  by  the  average  amount  of 
time  it  requires  for  service.  Thus  if  A^(t)  = N^(t  + dt)  - N^(t), 
i = 1,  2,  then  we  assume 

p N (t)/Vi 

P(A^(t)=0,  4j(t)-i|N(t))  = * °"5‘> 

The  motivation  for  using  these  transition  rates  is 
as  follows.  Over  a short  period  of  time,  say  (t,  t + dt) , at 
most  one  specific  job  has  positive  probability  of  being 
completed.  This  is  in  contrast  to  (3.2)  in  which  either 
type  has  positive  probability  of  completion.  Nevertheless,  if 
we  focus  on  a longer  period  of  time,  say  the  amount  of  time 
needed  to  service  all  jobs  currently  in  the  queue,  then  on 
the  average  N^{t)/iJ^  time  units  will  be  spent  servicing 
type  i jobs,  i = 1,2.  Thus  the  server  will  spend  approximately 
a fraction  (N^  (t ) /u^)  / (t) + N2(t)/vi2)  of  the  total 

time  servicing  type  i customers.  Approximating  (3.1)  by  (3.2) 
amounts  to  applying  this  fraction,  which  changes  dynamically. 
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to  every  interval  of  length  dt  rather  than  just  to  the 


longer  interval  in  which  all  customers  are  served.  This 
approximation  is  suggested  only  for  steady  state  calculations, 
since  it  does  not  model  the  local  behavior  of  a FCFS  server 
faithfully. 

The  transitions  (3.2)  belong  to  a general  parametric 
class  of  transitions  of  the  form 


UlN. (t) 

P(A^(t)=-l,  A2(t)=0|N(t))  = + cN2(t)  ^ 


(3.3) 


P(Aj^(t)=0,  A2(t)=-l|N(t)) 


VJ2cN  (t) 

(t)  + CN2  (t)  o(<3t) 


The  parameter  c governs  the  relative  weight  given 
to  jobs  of  either  type.  If  c = 1,  then  each  is  given  equal 
weight, which  corresponds  to  a processor  sharing.  If  c = Uj^/y2' 
then  the  queue  discipline  becomes  FCFS.  As  c -*•  +®  , type  2 
jobs  receive  priority,  while  if  c -►  0,  type  1 jobs  receive 
priority.  Equations  (3.3)  provide  a new  approach  to  assessing 
system  performance  as  a function  of  the  scheduling  or  queueing 
disciplines . 

Equations  (3.3)  coupled  with  the  usual  equations  govern- 
ing arrivals  to  the  queue  in  (3.1)  provide  the  transition  rates 
for  {(Nj^(t),  N2(t)),  t ^ 0}.  These  equations  can  be  solved 
numerically  using  Gauss-Seidel  iteration.  The  accuracy  of 
this  method  and  selected  numerical  results  will  be  given  later. 
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The  third  modelling  approach  utilizes  a diffusion 
approximation  based  on  (3.3).  This  method,  while  perhaps  less 
accurate  than  the  first  two  has  the  virtue  of  providing  simple 
closed-form  expressions  for  various  system  parameters  such  as 
the  expected  queue  lengths.  It  is  difficult  to  make  such 
assessments  using  numerical  methods. 

We  follow  the  method  outlined  in  Section  2 by  writing 
stochastic  differential  equations  to  approximate  (t) , 
i = 1,2.  Hence 

N,  (t) 

dNj^(t)  = Xp(N-N^(t)-N2(t))dt  - Vi^N  dt 


+ /Xp(N-Nj^(t)-N2(t)  ) dWj^(t) 


^1^  N]Tty+cNjTtT 


(3.4) 


N,  (t) 


dNjIt)  - X(l-p)(N-Hj(t)-Nj(t))dt  - MjN 


+ /X(l-p)  (N-Nj^(t)-N2(t)  ) dW3(t) 


4 


-\/UoN 


cr4"(t) 


2^  N^TtT+cNjTtT 


We  next  define  * (N^  (t) -Nq^  (t) ) //Jj,  i = 1,2 

so  that  (Nj^(t)  ,N2  (t)  ) = N (q^^  (t)  ,q2  (t ) ) + /N  (Xj^j^  (t)  ,X2j,  (t) ) . 
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The  first  term  provides  the  deterministic  approximation,  while 
the  second  is  the  stochastic  noise  term.  To  calculate  repre- 


sentations for  q^(t)  and  X^j^(t),  i = 1,2,  we  need  an 
asymptotic  expansion  for  (t) / (N^ (t) +N2 (t) ) , 


Nj^(t) 

N^TtJ+cNjTtT 


(t) 

q^Tty+cqJTty 


/N 


5 I + o ( — ) 

(qj^  (t)+cq2  (t)  ) / /N 


(3.5) 


cN^Ct) 

N^TtTT^NjTty 


cq2 (t) 

q^^  (t)  +cq2  (t) 


, c / ‘^l  ,l^ 

+ I ;; I + O ( ) 


/N  \ (qj^  (t) +cq2  (t)  ) 


/N 


We  now  let  N ■*■  + “ . Using  the  methods  outlined  in 
Section  2,  { (t)  ,X2j^  (t ) ) , t ^ 0}  converges  weakly  to 

{ (Xj^  (t)  ,X2  (t)  , t ^ 0)  provided  (qj^(t),  q2{t))  satisfies  a 
system  of  ordinary  differential  equations.  We  find  using  (3.4) 
and  (3.5)  that 


q{  (t) 


Xpd-q^^  (t)  -q2  (t)  ) 


Uj^q^^  (t) 
q^'(t)+cq2  (t) 


q^(t) 


X(l-p)  (l-qj^(t)-q2(t)) 


M2cq2 (t) 

q^TtT?cq^TtT 


(3.6) 


with  (qj^(O)  ,q2(0))  = (Nj^  (0)  ,N2  (0) ) /N. 

{ (Xj^  (t)  ,X2  (t) ) , t ^ 0}  will  satisfy 


The  noise  process 


dX(t)  = A.X.dt  + B.dW.  , 


X(0)  = 0 a.s. 


(3.7) 


where  X(t)  - (Xj^  (t ) ,X2  (t)  ) 


/-Xp  - 


MlCq^Ct) 


(qj^  (t)+cq2  (t)  ) 


\ 


-X(l-p)  + 


(qj^  (t)+cq2  (t) ) 


- Xp  + 


(t) 


(qj^  (t)+cq2  (t)  ) 


-X(l-p)  - 


(qj^  (t)+cq2  (t)  ) 


and 


/Ap(l-qj^(t)-q2(t)) 


?t=l 


/X  (l-p)  (l-qj^  (t)  -q2  (t) ) Vq^7ty+oqJ?tr 


Equation  (3.7)  defines  a bivariate,  nonstationary 
Ornstein-Uhlenbeck  process.  From  Arnold  (Chapter  8),  we  know 
E(X(t))  •=  0 for  all  t,  and  X(t)  will  have  a bivariate 

T 

normal  distribution  with  covariance  matrix  * E(X(t)  X (t)) 
given  by  the  unique  non-negative-definite  solution  of 
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0 


(3.8) 


it  “ \lt  ?t?t 


fo  = 


This  matrix  Riccati  equation  can  be  solved  numerically  (see 
Arnold,  Chapter  8 for  a general  expression) . 

The  mean  of  N(t)  is  N(qj^  (t)  ,q2  (t)  ) and  may  be  found 
by  solving  (3.6).  Because  (3.6)  is  non-linear,  it  will  have  to 
be  solved  numerically.  In  fact,  we  will  not  now  be  concerned 
with  the  dynamic  behavior  predicted  by  the  models  (3.6)  and 
(3.7).  Instead,  we  choose  to  examine  the  predicted  steady- 
state  behavior  that  is  obtained  by  letting  t -*^  + ®.  As 
t '*  snd  A^  A,  B.  We  find 

- i/(Pj^+  p^))/(i  + 


(3.9) 


qj  = (1  - l/(pj^  + P2)  ) p2/cPj_  (1  + P2/cpj^) 


where  - Xp/u^,  i = 1,2. 

N(t)  will  be  approximately  bivariate  normally  distributed 
with  mean  N(qj^,q2)  and  covariance  matrix  NI  where  E 
is  the  unique  non-negative  definite  solution  of 

AE  + Ea"^  = - Bb"^  (3.10) 


with 
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and 


(qi+cq2) 

(qi+cq2) 


Ap(l-qj-qj)  ♦ 


T 

BB^  = 


X(l-p) (l-q^-q2)  + 


One  use  of  (3.9)  and  (3.10)  is  to  study  the  influence  of 
the  parameter  c on  system  performance.  If  c = P2/Pj^»  then 
the  values  of  q^^  and  q2  are  equal.  If  c < P2/p2^f  then 
^2  ^ whereas  q2  ^ ^ p2/Pl*  should  be  noted 

that  no  matter  what  c is  chosen,  q^^  + q2  will  equal 
1 - l/(pj  + P2)  • This  is  intuitive,  since  q^^  + q2  represents 
the  number  of  jobs  in  service.  No  matter  which  value  of  c 
is  chosen,  the  same  backlog  of  jobs  will  be  faced  by  the  server. 
The  parameter  c controls  the  relative  fraction  of  jobs  in 
the  queue  of  each  type  and  thus  the  response  time  for  each 
job  type.  Here  “ ‘^Pi/p2*  costs  are  assigned  to 

waiting  or  turnaround  time,  then  an  appropriate  value  of  c 
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could  be  chosen  to  minimize  expected  cost.  This  value  of 
c can  then  be  implemented  by  using  the  corresponding  scheduling 
algorithms . 
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4.  Numerical  Results  for  Model  2 


In  this  section  we  present  numerical  results  to  assess 
the  accuracy  of  the  three  approximating  methods  given  in  Section 
The  methods  referred  to  in  the  following  tables  are: 

Method  1:  Exact  solution  of  (3.1)  using  Gauss-Seidel  iteration. 

* 


Method  2:  Calculations  made  using  u = + {l-p)/y2) 

and  assuming  the  number  of  Type  1 jobs  in  line  has 
a binomial  distribution. 

Method  3:  Weighted  processor-sharing  system  (3.2)  using  Gauss- 
Seidel  iteration. 

Method  4:  Diffusion  approximation  using  (3.2),  (3.9),  and  (3.10) 
Here  are  numerical  results  for  particular  cases  to  which  these 
methods  have  been  applied. 


3. 


Case  1.  N = 10, 

p = 1/2,  X = 10, 

= 

10, 

Ny2  = 20 

(.'. 

= 1/2,  P2  = 1/4) 

• 

Method  E(Nj^) 

E(N2) 

Var  (N^ 

) 

Var (N2) 

P(Nj^,N2) 

1 4.50 

4.17 

2.58 

2.41 

-.72 

2 4.50 

4.17 

2.47 

2.47 

-.73 

3 4.53 

4.14 

2.66 

2.29 

-.73 

4 4.33 

4.33 

2.70 

2.31 

-.73 

Here  p(Nj^,N2) 

refers  to  the  ordinary 

product-moment 

correlation 

between  Nj^  and  N2 . 
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10 


20 


Case  2.  N = 10,  p = 1/4,  A = 10,  = 

(•  • ~ r P 2 ~ ^/8)* 


, Nu2  = 


Method 

E(Nj^) 

E(N2) 

Var  (N^) 

Var (N^) 

P (Nj,N2) 

1 

2.25 

6.15 

1.79 

2.43 

-.62 

2 

2.25 

6.15 

1.73 

2.53 

-.63 

3 

2.28 

6.12 

1.84 

2.26 

-.61 

4 

2.10 

6.30 

2.07 

2.82 

-.69 

Case  3. 

N = 10,  p = 1/2,  A 

o 

II 

= 10, 

o 

in 

II 

CM 

(.*.  p^  = 1/2,  P2  = 

1/10) 

• 

Me  t hod 

E(N^) 

ECN^) 

Var(N^) 

Var(N2) 

P(Nj,N2) 

1 

4.50 

3.83 

2.67 

2.39 

-.56 

2 

4.50 

3.83 

2.39 

2.39 

-.64 

3 

4.58 

3.76 

2.83 

1.98 

-.63 

4 

4.17 

4.17 

2.99 

2.09 

-.66 

Case  4 . 

N = 10,  p = 3/4,  A 

= 10, 

= 10. 

Nu2  = 20 

(.*.  Pj^  = 3/4,  = 

1/8) 

• 

Method 

E(N^) 

E(N2) 

Var (N^) 

Var (N^ ) 

P(N^,N2) 

1 

6.75 

2.11 

2.33 

1.65 

-.71 

2 

6.75 

2.11 

2.24 

1.67 

-.73 

3 

6.76 

2.10 

2.40 

1.61 

-.73 

4 

6.64 

2.21 

1.98 

1.16 

-.63 
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These  four  cases  were  chosen  to  be  quite  extreme  in  order  to 
provide  a reasonable  test  of  the  various  methods.  We  have  also 
provided  the  variances  of  and  rather  than  the  standard 

deviations;  the  latter  would  exhibit  small  percentage  differences 
than  do  the  variances. 

For  N = 10,  there  would  seem  to  be  little  hope  that  Method  4, 
the  diffusion  approximation,  will  be  very  accurate.  This  method 
requires  N -►  + “> , yet  even  for  N = 10  the  accuracy  is  surprising 
especially  for  the  means.  In  all  cases  E(Nj^)  + E(N2)  is  very 
nearly  exact,  and  the  individual  components  are  off  by  at  most 
1/2.  This  error  arises  from  the  fact  that  the  different  types 
have  different  amounts  of  service  time  and  hence  different  server 
occupancies,  but  this  discrepancy  decreases  as  N increases. 

The  variances  are  less  satisfactory  in  terms  of  percentage  error, 
but  the  absolute  errors  are  probably  small  enough  to  be  practically 
negligible.  Accuracy  increases  with  N,  and  the  method  is 
appealing  because  of  the  closed  form  expressions  produced. 

Method  2 offers  a great  deal  of  accuracy  in  most  cases. 
Unfortunately,  it  is  very  special  and  is  not  easily  extended  to 
more  complex  models,  in  particular  to  those  with  a fixed  number 
of  terminals  of  each  type.  In  such  a case  one  must  keep  track  of 
the  order  of  jobs  in  the  queue.  Even  changing  from  an  exponential 
to  a general  phase  type  distribution  will  reduce  the  accuracy  of 
Method  2;  the  accuracy  of  variance  estimates  will  suffer  more 
than  that  of  mean  estimates. 
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* 


Method  3 gives  acceptable  results,  better  than  Method  4 but 


not  as  good  as  Method  2.  Nevertheless,  this  method  offers 
great  promise  in  that  it  can  be  generalized  to  handle  more  compli- 
cated systems  and  phase-type  service  distributions.  The  approxi- 
mation given  in  (3.3)  for  two  types  of  service  can  be  extended 
to  many  types  and  phase  type  service.  Furthermore,  it  offers 
the  way  to  studying  the  performance  of  the  system  as  a function 
of  the  queue  discipline.  If  N is  large  and  the  service  rates 
are  not  widely  different,  the  accuracy  of  Method  2 should  be 
accurate  enough  for  most  purposes. 
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5 . Model  3;  Two  Types  of  Terminals 


In  this  section,  we  remove  the  assumption  that  each  terminal 
can  submit  either  type  of  job,  a fraction  p of  which  are  of  Type  1. 
Instead,  we  assume  that  each  terminal  submits  only  jobs  of  a certain 
fixed  type.  Specifically,  let  terminals  submit  only  jobs  of 

type  i,  i = 1,2,  and  let  N = + M2  be  the  total  number  of 

terminals.  Jobs  of  type  i require  an  exponential  (y^)  process- 
ing time.  Later  we  will  generalize  to  allow  phase  type  service 
distributions . 

The  "fixed  terminal"  assumption  is  not  ideal  but  it  probabl/ 
represents  the  real  situation  more  faithfully  than  does  the 
assumption  that  terminals  can  submit  any  job  type  at  any  time. 

It  seems  plausible  that  terminal  users  typically  interact  with 
a single  problem  type,  hence  job  type,  for  a significant  period 
of  time.  As  a result,  the  fixed  terminal  assumption  is  likely  to 
be  reasonable,  if  not  ideal,  over  a moderate  period  of  time. 
Modifications  which  allow  for  the  sign-on  and  sign-off  of  users, 
changes  in  job  types,  and  the  influence  of  queue  discipline  on 
such  a system  will  be  considered  in  a subsequent  paper. 

We  let  N^(t)  denote  the  number  of  type  i jobs  in  service 
at  time  t.  The  arrival  rate  of  type  i jobs  to  the  computer  is 
A^(M^  - N^(t))dt  + o(dt),  i = 1.2. 

An  exact  analysis  of  this  stochastic  model  is  very  difficult 
because  the  order  of  customers  in  a FCF3  queue  must  be  kept  as 
a part  of  the  state  description.  The  pair  (Nj^  (t)  ,N2  (t ) ) alone 
is  not  sufficient  to  guarantee  the  Markov  property.  A Markovian 
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state  description  will  require  a very  large  state  space,  and  any 
solution  must  apparently  be  computed  numerically.  For  phase- type 
service  distributions  and  and  M2  around  25,  it  seems  that 

one  must  resort  to  simulation.  As  a result,  it  is  difficult  to 
assess  the  dependence  of  system  performance  variables  such  as 
queue  lengths  and  utilization  on  service  rates,  input  rates,  number 
of  terminals,  and  queue  discipline.  Consequently,  we  seek  approxi- 
mate solutions  which  lead  to  tractable  results.  We  present  three 
such  methods.  First,  the  weighted  processor  sharing  FCFS  approach 
introduced  in  Section  3 can  be  adopted.  This  will  not  lead  to 
closed  form  expressions  but  will  considerably  reduce  the  state 
space.  The  state  space  will  be  of  the  order  of  Mj^*M2  for 
exponential  service  times,  but  of  course  much  larger  for  phase- 
type  service  distributions.  Second,  we  give  a diffusion 
approximation  for  the  system  based  on  the  weighted  processor 
sharing  method.  Third,  a new  method  which  keeps  exact  track  of 
the  job  in  service  will  be  introduced.  The  third  method  has  the 
advantage  of  a small  state  space  even  for  general  phase  distri- 
butions and  is  more  accurate  than  the  other  two.  Even  so,  it 
is  designed  to  handle  only  FCFS  queues  and  does  not  seem  to  lead 
to  a simple  diffusion  approximation  with  associated  closed- 
form  expressions. 

The  weighted  processor  sharing  FCFS  approximation  is  defined 
by  the  following  transitions  and  transition  rates. 
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/ 


Transitions 


Rate 


(Nj  N^)  — > (N^-l.N^) 

(N^,N2-1) 

(Nj^  + 1,N2) 

(N^,N2+1) 


UiN^ 


Nj^  + cN^ 


dt  + o(dt) 


W2CN2 


+ CN2 


dt  + o(dt) 


(5.1) 


■*■  o(dt) 


X2(M2-N2)dt  + o(dt) 


One  can  derive  balance  equations  and  solve  them  numerically 


for  given  values  of  the  parameters  and  Uj 


Phase-type  distributions  can  be  easily  incorporated  by  defining 
the  number  of  type  i jobs  currently  in  phase  j of  their  service. 

. The  rate  at  which  such  jobs  complete  service  is  taken  to  be 
*^ij^ij*^ij  ^ constants  can  be  chosen 

to  model  various  types  of  queue  disciplines.  Generally,  letting 
c^j  ® give  FCFS,  c^^  * 1 will  give  processor  sharing, 

and  c^j  close  to  zero  or  infinity  represents  various  priority 
systems . 

To  obtain  more  insight  into  (5.1),  we  next  develop  a diffusion 
approximation  following  the  methods  outlined  in  Sections  2 and  3. 
Specifically,  we  treat  Nj^(t)  and  N2(t)  as  continuous  variables 


28 


whose  probabilistic  structure  is  given  by  the  stochastic 
differential  equations 


(t) 


Xl  (t)  )dt 


Nj^(t)  + cN2{t) 


dt 


+ /X^TMp-'N^Ttrr  dWj^(t) 


/ y,N  N (tl~ 
VNj^  (t)  + CN2  (t) 


(5.2) 


dN2(t) 


M-N  cN-  (t) 

^2  (M2-N2  (t)  )dt  - Nj^(t)  + cN2(t) 


+ /X^TMj^^j^TETT  dW^Ct)  - 


~M^N~cNjTtl 

Nj^(t)  + cN2(t) 


dW^ (t) 


We  next  define  X^j^(t)  = (N^(t)  - Nx^(t))/»^,  i = 1,2  and 
let  = Nj^(t)/N,  i = 1,2.  Again,  if  (Xj^(t),  X2  (t)  ) satisfies 

a certain  set  of  ordinary  differential  equations,  then  as  N -►  <» 

t ^ 0}  converges  weakly  to  a limiting  stochastic 
process  {{Xj^(t),  X2(t)),  t ^ 0}.  The  deterministic  approximation 
is  given  by  N(Xj^(t),  X2(t))  where 

x|(t)  = (t)  ) - Uj^Xj^  (t)/(Xj^  (t)  + cx2(t)) 

(5.3) 

x^(t)  - A2  (oi2”*2  ^ ■ ^2^X2  (t)/(Xj^  (t)  + cX2(t)) 
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«f 

I 


with  Xj^(O)  * Nj^(0)/N.  The  noise  process  {(Xj^(t),  X2  (t) ) , t ^ 0} 
will,  according  to  Ito's  Lemma,  satisfy 

dX(t)  = A^X(t)dt  + B^dW(t)  , X(0)  = 0 (5.4) 

where  X(t)  = (Xj^  (t)  , X2(t))'^  , 


The  t-dependence  of  and  X2  has  been  suppressed  in 

and 

Equation  (5.5)  characterizes  a nonstationary  bivariate 
Ornstein-Uhlenbeck  process  with  mean  0.  The  covariance  matrix 
of  (Nj^  (t ) ,N2  (t) ) is  thus  approximately  given  by  where 

is  the  unique  non-negative  definite  solution  of  (3.8).  The 
mean  of  (N^^  (t ) ,N2  (t ) ) is  approximately  N (Xj^  (t)  , X2  (t) ) . 
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The  steady-state  behavior  of  the  system  is  deduced  by 


letting  t -►  <».  The  steady  state  means,  (Xj^,X2),  will  satisfy 


0 = - Uj^Xj^/CXj^  + CX2) 

0 = X2(°‘2”^2^  “ ^^2^  ' 


(5.5) 


while  the  covariance  matrix,  NE^,  can  be  deri /ed  from  (3.10). 

This  analysis  is  presented  only  for  exponential  service, 
but  can  be  easily  extended  to  allow  for  general  phase-type 
service  distributions.  Numerical  examples  will  be  presented  in 
the  next  section. 

The  third  approximation  method  offers  excellent  accuracy. 

In  this  method  the  state  description  keeps  track  of  the  number 
of  jobs  of  each  type  awaiting  service,  as  well  as  the  phase  of 
the  current  job  receiving  service.  The  exact  order  of  jobs 
awaiting  service  is  not  recorded.  Instead  when  a customer 
finishes  service,  the  next  customer  is  chosen  at  random  from  those 
waiting  for  service.  For  exponential  service,  the  state  of 
the  system  is  given  by  {(Nj^(t),  N2(t),  I(t)),  t ^ 0}  where 
N^(t)  represents  the  number  of  type  i jobs  awaiting  service 
and  I(t)  gives  the  identity  of  the  job  in  service.  Then  the 
transitions  are  given  by 
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Transition 


Rate 


(Nj^,N2,1) 


(Nj^,N2,2) 


— > (Nj^+1,N2,1) 
(Nj,N2+1,1) 

(N^-1,N2,1) 

(Nj^,N2-1,2) 


— > (Nj^+1,N2,2) 

(Nj^,N2+1,2) 
(Nj^-1,N2,1) 
(Nj^,N2-1,2) 


Xi(Mi-Nj^)dt  + o(dt) 
Xj^(M2-N2)dt  + o(dt) 
li^N(Nj^/(Nj^+N2)  )dt  + o(dt) 
jjj^N(N2/(Nj^+N2))dt  + o(dt) 


Xi(Mi-Nj^)dt  + o(dt) 
X2(M2-N2)dt  4-  o{dt) 
y2N(Nj^/(Nj^+N2)  )dt  + o(dt) 
M2N(N2/(Nj^+N2)  )dt  + o(dt) 


(5.6) 


This  approximation  can  be  easily  extended  to  allow  for 
phase-type  service  distributions,  but  apparently  does  not  lend 


itself  to  closed-form  approximations. 
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6. 


Numerical  Results  for  Model  3. 


In  this  section,  we  present  numerical  results  for  assessing  the 
accuracy  of  the  three  approximation  methods  given  in  Section  5.  The 
methods  referred  to  in  the  following  tables  are: 

Method  1:  Solution  of  the  exact  system  using  simulation. 

Method  2;  Approximation  based  on  (5.8)  using  simulation. 

Method  3:  Approximation  based  on  (5 , 1 ) -weighted  processor  sharing 
with  c = VJj^/P2  using  simulation. 

Method  4:  Diffusion  approximation  using  (5.5)  and  (3.10). 


For  the  first  three  methods,  a simulation  consisting  of 
600,000  system  transitions  blocked  into  20  groups  of  30,000  was 
performed.  The  estimated  standard  deviation  is  given  in  parentheses. 


All  cases  consisted 

of  Gamma 

(2, 

Np^)  service 

times 

with 

>^1  = 

1. 

. , = 1/2 

and  Mj^  = 

M2  = 20. 

Case  1. 

X — 2 , X 2 

= 2 

Method 

E(Nj^) 

E(N2) 

Var (N^) 

Var (N^) 

p(Nj^ 

1 

16.64  (.05) 

16.67  ( 

.05) 

3.17 

(.15) 

2.46 

(.07) 

.02 

(.03) 

2 

16.65  (.07) 

16.65  ( 

.04) 

3.19 

(.10) 

2.47 

(.06) 

.02 

(.02) 

3 

16.66  (.05) 

16.67  ( 

.09) 

3.20 

(.14) 

3.20 

(.14) 

.12 

(.03) 

4 

16.67 

16.67 

2.97 

3.14 

.09 
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Case  2 


1 


Case  2. 

* 

2,  X2 

= 1 

Method 

E(N 

ll 

E(N2) 

1 

16.23 

(.06) 

13.72  (.10) 

2 

16.25 

(.08) 

13.72  (.10) 

3 

16.26 

(.06) 

13.74  (.11) 

4 

16.28 

13.72 

Case  3. 

^1  = 

1,  X2 

= 2 

Method 

E(N 

ll 

E(N2) 

1 

14.01 

(.08) 

16.49  (.05) 

2 

14.02 

(.10) 

16.48  (.04) 

3 

14.03 

(.09) 

16.49  (.06) 

4 

14.03 

16.49 

Case  4 . 

^1  = 

1,  X2 

= 1 

Method 

E(N 

E(N2) 

1 

13.31 

(.08) 

13.34  (.10) 

2 

13.28 

(.12) 

13.33  (.09) 

3 

13.31 

(.10) 

13.36  (.11) 

4 

13.33 

13.33 

Var(Nj^) 

Var(N2) 

p(Nj^, 

N2) 

3.53 

(.13) 

4.22  (.15) 

.06  ( 

.03) 

3.54 

(.18) 

4.15  (.17) 

.07  ( 

.02) 

3.46 

(.11) 

5.41  (.04) 

.19  ( 

.04) 

3.36 

5.45 

.16 

Var  (Nj^) 

Var (N2) 

P(N, 

X 

.Nj) 

4.77 

(.18) 

2.57 

(.08) 

.06 

(.02) 

4.91 

(.23) 

2.55 

(.09) 

.05 

(.02) 

4.75 

(.15) 

3.37 

(.12) 

.14 

(.03) 

4.69 

3.33 

.11 

Var  (Nj^) 

Var(N2) 

P (N^ 

,N2) 

5.35 

(.18) 

5.41 

(.22) 

.12 

(.03) 

5.53 

(.27) 

4.46 

(.24) 

.11 

(.03) 

5.38 

(.19) 

5.87 

(.30) 

.23 

(.03) 

5.28 

5.88 

.20 
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Case  5. 

= 

1,  ^2 

= 1/2 

Method 

E(N 

li 

E(N2) 

1 

11.60 

(.12) 

8.33 

(.12) 

2 

11.63 

(.12) 

8.35 

(.14) 

3 

11.62 

(.14) 

8.41 

(.19) 

4 

11.72 

8.28 

Case  6 . 

^1  = 

1/2, 

X2  = 1 

Method 

E(Nj^) 

E (N 

2^ 

1 

9.27 

(.10) 

12.68 

( .11) 

2 

9.28 

(.13) 

12.68 

(.10) 

3 

9.26 

(.14) 

12.70 

(.12) 

4 

9.28 

12.68 

Case  7. 

^1  = 

1/2, 

^2  = 1/2 

Method 

E(Nj^) 

E(N2) 

1 

6.60 

(.18) 

6.74 

(.14) 

2 

6.62 

(.21) 

6.75 

(.20) 

3 

6.59 

(.16) 

6.82 

( .17) 

4 

6.67 

6.67 

Var  (Nj^) 

Var (N2) 

p(Nj^ 

,N2) 

6.74 

(.30) 

6.13 

(.35) 

.26 

(.03) 

6.87 

(.40) 

6.05 

(.35) 

.27 

(.03) 

7.37 

(.31) 

7.64 

(.55) 

.40 

(.04) 

7.02 

8.12 

.37 

Var  (Nj^) 

Var(N2) 

P(N^ 

,N2) 

6.27 

(.22) 

5.04 

(.23) 

.20 

(.02) 

6.37 

(.31) 

4.90 

(.22) 

.19 

(.02) 

6.46 

(.31) 

6.67 

(.33) 

.27 

(.03) 

6.36 

6.64 

.25 

Var(N^) 

Var(N2) 

P(N^,N2) 

7.57 

(.33) 

7.11 

(.31) 

.42  (.03) 

7.67 

(.57) 

7.11 

(.50) 

.43  (.04) 

8.29 

(.43) 

9.00 

(.54) 

.52  (.03) 

8.40 

9.88 

.52 
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Case  8.  = 1/4,  = 1/2 


Method  E(Nj^) 

1 2.71  (.11) 

2 2.69  (.14) 

3 2.76  (.12) 

4 2.46 


E(N2) 

4.86  (.17) 
4.81  (.19) 
5.09  (.15) 
4.38 


Var (N^) 

4.05  (.22) 
4.12  (.22) 

4.28  (.34) 

5.28 


Var (N2) 

7.31  (.37) 
7.31  (.32) 
9.12  (.34) 
13.04 


f (Nj^,N2) 

.49  (.02) 
.50  (.02) 
.54  (.02) 
.66 


The  numerical  examples,  which  cover  a wide  range  of  traffic 
intensities,  illustrate  the  extraordinary  agreement  of  the  three 
methods  in  estimating  the  means.  Even  the  diffusion  approximation 
offers  nearly  exact  answers.  For  the  variances.  Method  2,  which 
keeps  track  of  the  customer  type  in  service,  appears  to  be  nearly 
exact.  The  differences  between  1 and  2 are  dominated  by  the 
quoted  error  of  simulation.  On  the  other  hand.  Methods  3 and  4 
do  not  provide  good  accuracy,  especially  for  Var(N2)  and  p(Nj^,N2). 
It  is  perhaps  surprising  that  the  diffusion  approximation  (Method  4) 
agrees  so  well  with  Method  3 in  view  of  the  small  value  of  N 
(N  ■ 40).  It  is  clear  that  for  phase-type  distributions,  the 
variances  reported  by  the  weighted  processor  sharing  FCFS  approxi- 
mation will  be  inaccurate,  providing  only  a rough  estimate  of 
the  true  value.  On  the  other  hand,  the  means  are  exceptionally 
good,  and  the  diffusion  approximation  means  provide  trustworthy 
expressions  for  average  system  occupancy  of  each  type. 
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We  have  not  shown  here  the  computer  idleness  probability. 
It  behaves  similarly: 

Method  2 provides  nearly  an  exact  result  while 
Method  3 gives  about  10  to  20  percent  error. 

The  diffusion  approxiir.=»tion  applied  here  does  not  recognize 
idleness  of  the  server. 
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